Constraining density functional approximations to yield self-interaction free potentials 
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Self-interactions (Sis) are a major problem in density functional approximations and the source of 
serious divergence from experimental results. Here, we propose to optimize density functional total 
energies in terms of the effective local potential, under constraints for the effective potential that 
guarantee it is free from SI errors and consequently asymptotically correct. More specifically, we 
constrain the Hartree, exchange and correlation potential to be the electrostatic potential of a non- 
negative effective repulsive density of A'^ — 1 electrons. In this way, the optimal effective potentials 
exhibit the correct asymptotic decay, resulting in significantly improved one-electron properties. 



I. INTRODUCTION 

It is already thirty years since Perdew and Zunger 
in a seminal contribution [l| proposed to cure the self- 
interaction (SI) error in density functional approxima- 
tions (DFAs). The SI error arises from the incomplete 
cancellation of the self-repulsion of the electron density p 
in the direct Coulomb or Hartree energy U\p\ by the ap- 
proximate exchange energy functional E^ [p]. SI errors 
manifest in inaccuracies of DFAs [2| in many ways, e.g. 
in the calculation of binding energies [3J, underestimation 
of activation energy barriers [j, |5|, and in single-particle 
properties like ionization potentials (IPs) [a, |7| , electron 
affinities (EAs) (unbound anions) Q and band gaps of 
solids @- 

Perdew and Zunger [l| proposed a many-body SI cor- 
rection energy term which, in the limit of a single elec- 
tron, eliminates the SI error exactly. Their work initiated 
the field known as self-interaction corrected density func- 
tional theory (SIC-DFT). Unfortunately, the many-body 
generalization of the one-electron SI energy correction is 
not unique and to date an unambiguous definition is not 
available. Rigorously, we have a sufficient condition for 
an approximate exchange and correlation (XC) energy 
density functional Ey^c [p] to be A^-representable [10| and 
thus free from many-body SI errors. A distinction of 
one- and many-bod y S I error was made independently 
by Ruzsinszky et al[ll|. An important development in 
this area is the appreciation of the underlying relation- 
ship between the SI error with the fractional charge error 
[l|, |6|, llO| . The approximate treatment of many-body SI 
errors with SIC-DFT leads traditionally to single-particle 
equations with orbital dependent potentials, i.e. the min- 
imization problem is significantly more complicated than 
an iterative diagonalization. For solids, SIC-DFT is ex- 
pressed in terms of maximally localized Wannier states. 
An advantage of removing SI errors is that it improves 
orbital energies p^ll3l|. These orbital energies, are com- 
monly obtained as the eigenvalues of the non-diagonal 
Lagrange multiplier matrix employed to enforce orbital 
orthogonality, although the diagonal values of this ma- 
trix have also been proposed as appropriate [ij] . Despite 



complications, SIC-DFT has been extensively developed 
and applied to a large variety of systems pTl Il5l - [l9l | . 



A constrained minimization of the total energy, so that 
the occupied Kohn-Sham (KS) orbitals arise as the lowest 
cigcnorbitals of a common local multiplicative potential 
has also been developed recently [20|[, usiiig the optimized 
effective potential (OEP) methodjl, [2ll-[p] . 



Probably the most serious flaw caused by SI errors lies 
in the asymptotic behavior of the KS potential [2^ . If 
the cancellation of SI terms was complete then at infin- 
ity, the electron-electron contribution to the KS potential 
(Hartree and XC) should be {N — l)/r where N is the 
number of electrons. The physical meaning is obvious, 
at infinity each electron feels the screening of the nuclear 
charge by the remaining A'^— 1 electrons. The components 
of the Hartree potential ■^h and of the exact XC potential 
fxc to the asymptotic decay are N/r, and — l/r, respec- 
tively. However, the asymptotic decay of fxc"^ in typical 
DFAs, such as the local density approximation (LDA) 
or the generalized gradient approximation (GGA), does 
not follow a power law (— c/r), but is exponentially fast 
(c = 0). Consequently an electron at infinity is repelled 
by an effective charge of N rather than A — 1 electrons. 
The incorrect asymptotic behavior has dramatic conse- 
quences on one-electron properties like the IPs, EAs and 
the fundamental gaps. It also impairs the optical spec- 
trum through linear response in time dependent DFT|25| . 



In the present work, we propose a quantification of the 
many-body SI error in the Hartrec-exchange and corre- 
lation potential (fnxc) of any DFA (Sec. |Tl|. Following 
our definition, wc develop a simple way to address Sis 
by a constrained OEP minimization of the DFA total 
energy (Sec. llll| ). The constraints on the optimized effec- 
tive potential remove the effects of Sis from the potential 
and enforce the correct asymptotic behavior. Finally, in 
Sec. IIV| we present our implementation and the first nu- 
merical results. 



II. QUANTIFYING SELF-INTERACTIONS IN 
THE EFFECTIVE POTENTIAL 



Aiming to address in an unambiguous manner the Sis 
in finite systems we decided, rather than deahng with 
the approximate Hartree (U) and Exc energies which 
remain unchanged and contaminated with SI errors, to 
focus on the cfFcctivc local potential. The latter (i.e. the 
HXC potential) is obtained as the functional derivative 
of the HXC energy with respect to the density. In KS 
theory, iihxc = ^h + I'xc screens the nuclear attraction 
felt by a KS electron. By virtue of Poisson's equation, 
and following Gorling [23l , the Laplacian of whxc defines 
the HXC density phxc- 



V^WHXc(r) = -47rpHXc(r) 



(1) 



Here, pnxc is the density with electrostatic potential 
whxc- The integrated charge. 



Qhxc = dr PHXC (r) , 



(2) 



allows us to quantify the SI error of the approximate 
HXC potential. For example, if Qhxc = N —1 then each 
electron interacts with an effective electrostatic charge of 
A'' — 1 electrons and this is a necessary condition that the 
approximation is SI free. If A^ — 1 < Qhxc < N there is 
partial cancellation of Si's and finally if Qhxc = N, each 
electron interacts electrostatically with the charge of the 
other iV — 1 electrons plus an additional electron that 
can only be attributed to the same electron itself. We 
say that the corresponding HXC potential exhibits full 
SI effects. In popular approximations, such as LDA or 
GGA, Qhxc = A^ as can be deduced from the asymptotic 
behavior of the potential. Thus, according to the present 
definition there is full SI in the HXC potential in these 
approximations. 

Our criterion for SI can be equivalently expressed in 
terms of the exchange and correlation charge, 

Qxc = Ohxc - N , 

which must equal Qxc = ^1 for SI- free approximations. 
This condition was employed by Gorling j27| to constrain 
the asymptotic behavior of a finite-basis set implementa- 
tion of the exact exchange (EXX) potential. 

The exchange and correlation charge Qxc and the 
corresponding exchange and correlation density [23, [S^ , 
Pxc(r) = PHXc(r) — p(r), bear a similar name to the fa- 
miliar exchange and correlation hole of an electron at r, 
nxc(r, r')[2g|. The latter is a property directly obtained 
from the pair correlation function h{r, r'), 



nxc{r,r') ^ p{r') h{r,r') , 



(3) 



and satisfies the sum rule /dr' nxc(r,r') = —1. 

There is no easy and direct relation between the ex- 
change and correlation hole (two-electron property) and 



the exchange and correlation density (one-electron prop- 
erty). In particular, the approximate LDA XC hole. 



-.LDA 



(r, r'), satisfies the sum rule / dr' 



'■XC 



,r') = -l, 



since LDA corresponds to a physical system, the uniform 
electron gas, but the satisfaction of this sum rule does 
not preclude SI errors from the LDA potential, vl^^, ex- 
cept of course when LDA is applied (tautologically) to 
the uniform electron gas [26| . 



III. CONSTRAINING THE OPTIMIZED 
EFFECTIVE POTENTIAL 



In order to correct Sis in DFAs we propose to replace 
the Hartree, exchange and correlation potential f^xc ^^ 
the KS equations with the effective repulsive potential, 
Wrcp, obtained from a constrained minimization of the 



DFA total energy, 



Ucn(r) +Urop(r) 



0,(r) ==e,0,(r), (4) 



where Won(r) is the attractive electron- nuclear potential. 
The potential Wrop is optimized, in the fashion of the 
OFF method, by requiring that its N lowest orbitals give 
the density p = J2i=i \4>i\'^ that minimizes the DFA total 
energy. We constrain the DFA total energy minimization 
by restricting the potential Uiop to satisfy two conditions 
on its effective repulsive density Prcp, i.e. the density 
with electrostatic potential Wi-op(r), 



t^rcp(r) = / dr 



I Prcp(r') 



The two conditions read: 



?rcp ^ dr Prcp(r) = iV - 1 



Pr. 



> 0. 



(5) 

(6) 
(7) 



The system of equations ([6]) and ([7]) is equivalent to 
the numerically simpler-to-implcment system (jS]) and ([5|) , 
with 



f dr \p,,p{r)\ = N ■ 



1. 



(8) 



We remark that our constrained minimization results in 
a DFA total energy minimum that is in general higher 
than the global minimum, unless the global minimizing 
potential happens to satisfy the two extra constraints (|6]) 
and dZl). 

As wc have discussed, the normalization of the total 
effective repulsive charge ([6]) is necessary for the absence 
of Sis from the effective repulsive potential. However, 
Eq. © on its own is not sufficient to ensure the absence 
of SI effects in Vycp, as it would be energetically favorable 
to retain Sis locally near the system, plus a compensat- 
ing negative charge far away from the system to satisfy 
(pl). Specifically, by imposing this constraint alone, one 



will still obtain, almost everywhere, the global minimum, 
i.e., Vrep = i'hxC' ■^ith Prep = Phxc almost everywhere, 
integrating to almost TV up to a very large distance from 
the system, plus a compensating electronic charge of —1 
distributed over a large radius, very far away from the 
system, where the addition of the extra negative charge 
will not cost energetically. Since the large radius must 
be as far away as possible and its position is not well 
defined; the problem strictly has no solution. 

With the second (strong) condition, Eq. ([7]) or ([5]), 
such a pathological behavior can be avoided and the pair 
of conditions (|6]) and ([7]) becomes sufficient for the ab- 
sence of Sis from the potential - although probably it is 
not necessary any more. Equation ([7]) is an approximate 
condition but it has an obvious and appealing physical 
interpretation: since prep is non-negative and integrates 
to A^ — 1, it corresponds to a virtual system of A^ — 1 
electrons repelling the electron at r. There is no longer 
any freedom for our solution to collapse to the global 
minimum solution with a compensating negative charge 
at large distances. Hence the constraint of Eq. ([7]) allows 
for a physical enforcement of the constraint of Eq. © 
correcting Sis from Uj-cp- 

We note that if the effective repulsive density is ex- 
panded in a small finite basis set it is possible that em- 
ploying (|6]) alone may result in a solution that appears 
physical. This is an artifact of the smallness of the basis 
set, and by increasing the size of the effective repulsive 
density basis set, the pathology of not having a sufficient 
condition will emerge. 

In contrast to ug^^ ^ (5(t/[p]+£'xc^[p])/(5p, the poten- 
tial Wrop is not the functional derivative of t/ [p] -h £'xc'^ [p] 
and it is unknown if Wrcp is the functional derivative of 
some other functional. This question will be explored in a 
future work where our method is generalized to extended 



systems. Of course, the full KS potential Vg 



^rep 



in Eq. 



((4]) is the functional derivative with respect to the den- 
sity of the non-interacting kinetic energy functional Tg [p] 

8^^ P — J2i=i l0iPj as with any OEP theory. 

In a related work to correct the asymptotic behavior 
of the effective potential, Wu et al [2^ partitioned the ef- 
fective potential into the Fermi Amaldi potential, which 
has the correct asymptotic behavior, and they expanded 
the remainder in a finite localized basis set. However, for 
a large basis set (not complete), the tail of the potential 
(for moderate distances) will revert to that of the uncon- 
strained DFA potential and only for very large distances 
will the correct asymptotic behavior be recovered. 

Andrade and Aspuru-Guzik [30|, also aimed at the 
same problem using an ad hoc correction of the XC den- 
sity at large distances. In our approach the effective re- 
pulsive density and potential are obtained directly from 
the minimization procedure. 

A feature of our method which contrasts it to SIC- 
DFT, is that the DFA total energy is unchanged and con- 
sequently it is invariant under unitary transformations of 
the occupied orbitals. 

Finally, an important test for any theory that corrects 



Sis is that it has the correct one-electron limit. In this 
case, the repulsive potential Wiep should vanish for a one- 
electron system and indeed it is straightforward to con- 
firm that constraints (O and ([7]) give /9rcp(r) = 0, leading 
to Wrop(r) = 0, as expected. 

Our search for the effective repulsive density and po- 
tential is performed by expanding them in a basis and 
then searching for the expansion coefficients. 



.(r) 



/ 



viXi(r) 



(9) 



«rcp(r) ^Y."^Ur), with CKr) = / rfr' ^^ (10) 



where Xi(r) is an auxiliary basis set, for example, local- 
ized gaussians. 

The minimization of the approximate total energy, 
-Ed FA, with respect to vi under the conditions (|6]) and 
([8]) leads to the variation equation 



dE, 



DFA 



dvi 



fiXi + XXi, 



(11) 



where ^, A are Lagrange multipliers to satisfy (|6]), (|8]) 
and 



^1=1 d.rxi{r), 
Xi ^ dr xiir) 



I Prep (r) I 
Prcp(r) 



(12) 
(13) 



The derivative on the left hand side is obtained, as in 
the OEP method, through a chain rule: 



dE, 



DFA 



dvi 



drdr'p^P^l^^^. (14) 
6v,.cp{r) dprcp(r') dvi 



The functional derivative (Si?DFA/'5wrcp(r) is obtained 
analogously to the OEP functional derivative. 



SE^ 



DFA 



<^Wrop(r) 



-2E- 



^>Hxc ^^'^p\-^ ^^^r)Mr) , (15) 



where Whxc ^ ^i^ip] + ^xc^ip\)/^P' * ^^^s over occu- 
pied, a over unoccupied eigenorbitals of Urop, with e^ and 
Sa the corresponding eigen-energies. We also have that 



1 



^Prop(r') |r-r'| 



(16) 



We finally obtain 
SEdfa 



Svi 



HXC _ rep 



la^ (17) 



Sl'j = I drUr)Mr)^iir). 



(18) 



Here, vl^^ and v^'^ are the matrix elements of the po- 
tentials Wiep and WhxC' respectively. Eqs. pT|) and pT)) 
define a non-linear system of equations with respect to 
vi. This system can be solved, using an iterative scheme 
of two steps. In the first step, a linear system is solved 

by keeping the quantities (^i, 0q, e^, ea, vf^'~^, and Sj 
frozen. In the second step, a single-electron Hamiltonian 
problem is solved with the potential obtained in the pre- 



vious step and the quantities 



„HXC 



^15 V^a: ^z; ^ai ^ia 



and 



SlJ are updated. Our numerical implementation proved 
that this scheme is very efficient and usually only a few 
iterations are required to converge using a mixing scheme 
similar to Kohn-Sham iterative procedure. The potential 
obtained at the first step, i.e. when orbitals are frozen to 
the Kohn Sham orbitals is already a very good approx- 
imation to the local potential. The linear system that 
needs to be solved in each iteration has the form 



/ 



Akivi = bk+ fJ.Xk + XXk, 



with 

Aki 



E 






and bki = 



c(fc) ,,HXC 



(19) 



(20) 



The Lagrange multipliers /i, A are given by the solution 
of a simple 2x2 linear system obtained from Eq. (fTQ]) by 
multiplying both sides by the inverse of A, then by X^ 
(or Xk) and summing over fc, and using Eqs. ([5] and [5]). 
Eqs. (|19p and ((20|) constitute a simple modification of 
the usual OEP equations. The solution of Eq. (fT9|) is 
complicated because often the matrix A is singular. This 
problem is well known in the OEP method with finite 
basis sets J3l| - l35| , and the solution involves the inversion 
of A in the space of its non-singular eigenvectors, usu- 
ally with a singular value decomposition (SVD). How- 
ever, even after the SVD, and depending on the partic- 
ular basis sets, the resulting effective potential may look 
unphysical. In Ref. l36l we argue that in addition to the 
known technical problem of inversion of A lies an unex- 
pected discontinuity of the optimal potential, when the 
orbital basis set is truncated with a finite basis. In the 
present work, the effect of this discontinuity is reduced 
significantly by the restriction of the admissible poten- 
tials to satisfy conditions © and especially ([S]). 



IV. NUMERICAL APPLICATIONS 

Our numerical implementation is based on a Gaussian 
basis set expansion for the orbitals as well as for the 
effective potentials or for the effective repulsive density. 
XC functionals are provided by the Libxc librarvjSq. 



exact exchange functional and compare the potentials ob- 
tained with and without the constraints of Eqs. ^ and 
(O. This comparison is shown in Fig. [T]for the Ne atom. 
To obtain the unconstrained potentials in the two plots, 
the Hartree and exchange potential (i.e. not the effective 
repulsive density) was expanded directly in the uncon- 
tracted cc-PVTZ and uncontracted cc-PVQZ basis sets 
respectively. For the constrained case, the effective repul- 
sive densities were expanded in the same basis sets. The 
orbitals were expanded in the cc-PVTZ and cc-PVQZ 
basis sets. 

In the unconstrained minimization case, to obtain rea- 
sonable EXX potentials with the finite basis sets, we em- 
ployed the amended finite basis OEP Eq. (40) in Ref. [Sg 
which contains an extra term (discontinuity correction, 
with A = 10"'^) that restores continuity of the poten- 
tial. The unconstrained potentials contain an arbitrary 
constant and were shifted so that the energy of the high- 
est occupied molecular orbital (HOMO) equals that of 
Hartree Fock (HF). As can be seen in Fig.[l] the uncon- 
strained finite basis EXX potentials are very close to the 
full numerical EXX potential from Ref. [33l 

For the constrained minimization case, the discontinu- 
ity of the OEP with finite basis sets is reduced because 
the variational freedom of the admissible potentials is re- 
stricted by the two constraints. Then, inclusion of the 
discontinuity correction of Ref. [3^ is not necessary to ob- 
tain smooth potentials. However, subtle features of the 
EXX potential, such as the shell structure do not appear 
without the more complete description of OEP including 
the discontinuity correction. Thus, for a meaningful com- 
parison between the two results we show the constrained 
EXX potentials employing the discontinuity correction 
with A = 10~^. The two constrained potentials are al- 
most on top of the unconstrained potentials and the full 
numerical solution as seen in Fig. [T] In addition, the IP 
of Ne, given as the minus of the energy of the HOMO of 
the constrained potential, is almost identical to that of 
HF theory (^^0.2 eV lower) with the non-local exchange 
potential. This should be contrasted to the dramatic ef- 
fect of our constraints on the IPs in the case of LDA, i.e. 
a DFA with full Sis, as we will see in Sec. IIVBI 

To summarize this test, our results indicate that the 
constraints introduced here do not have a significant ef- 
fect on the potential for theories that are free from Sis. 
In Sec IIVBI we apply our methodology to LDA as an 
example of a DFA that is contaminated with SI errors. 
As we shall see in Fig. ^ the extra constraints modify 
the LDA KS potential substantially. 



B. Constraining the LDA potential 



A. Constraining the EXX potential. 

A rigorous test is to apply our constrained OEP 
method to a functional that is free from Sis such as the 



To illustrate our approach we chose the LDA 
functional [37| and we refer to the combined method as 
constrained LDA (CLDA). The LDA-KS potential misses 
the shell structure present in EXX even after the imple- 
mentation of our SI correction. For CLDA, the slight 
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Figure 1: Exchange part, ux, of the effective repulsive potential for the Ne atom using Eqs. ^ and ((Tjl (constrained potential) 
compared to the unconstrained finite basis EXX potential and the full numerical result or Ref. l3a . For the constrained and 
unconstrained EXX potentials, we used our amended finite basis OEP equations in Ref. ISb . 
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Figure 2: The XC part, uxc, of the CLDA-effective and 
KS-LDA potentials as well as the exchange part of the full 
numerical EXX potential (from Ref. [3j) for the Ne atom. 
— 1/r is also shown. In the inset, the effective repulsive charge 
density is shown. The notation (X-Y) stands for cc-pVXZ and 
uncontracted cc-pVYZ for the expansion of the orbitals and 
the effective repulsive charge density respectively. 



improvement in accuracy afforded by the inclusion of the 
discontinuity correction of Ref. |3a was not sufficient to 
warrant its use. We found that a simple SVD and inver- 
sion of the matrix A in Ec^. p9p gave adequate accuracy. 

To perform the SVD of A, we divide the space spanned 
by the eigenvectors of A in the null space and the rest, 
using a small parameter as a cut-off for the null eigen- 
values. For the systems and basis sets we considered, 
6 = 10^^ has proven a reasonable choice. 

In Fig. [21 we show the LDA and CLDA potentials for 
the Ne atom using finite basis sets as well as the effective 
repulsive density. Evidently, the effective potential ob- 
tained with CLDA has the correct asymptotic behavior. 
Also, the effective repulsive density is converged with re- 
spect to the different basis sets. Potentials with the cor- 





Basis 


Qu 


AE 


IP (LDA) 


IP(CLDA) 


Exp 


He 
Be 

Ne 

H2O 

NH3 

CH4 

C2H2 

C2H4 

CO 

NaCl 


T-Q 
T-T 
T-T 
T-T 
T-T 
D-D 
D-D 
D-D 
D-D 
D-D 


3.0 

6.0 
6.0 
1.5 
1.9 
3.9 
2.0 
1.2 




10-'* 


10-^ 
IQ-'^ 
10"^ 
10-^ 
10-^ 
lO"'^ 
10-2 


1.5 
1.1 
2.7 
1.1 
8.2 
2.7 
4.1 
1.1 
3.6 
6.8 


10"^ 
10-* 

lo-'' 

10-^ 
10-'' 
10-* 
10-^ 
10-^ 
10-* 
10-* 


15.46 
5.59 
13.16 
6.96 
6.00 
9.28 
7.02 
6.67 
8.75 
5.13 


23.14 
8.62 
18.94 
11.24 
9.81 
12.52 
10.63 
9.57 
12.73 
7.87 


24.6 
9.32 
21.6 
12.8 
10.8 
14.4 
11.5 
10.7 
14.1 
8.93 


F" 

cr 

OH- 

NH2" 
CN- 


rpa rp 
rpa rp 
rpa rp 
rpa rp 
rpa rp 


1.0 
1.0 
4.0 
4.0 
1.0 


10-* 
10-^ 
10-^ 
lO"'^ 
IQ-'^ 


2.7 
1.6 
1.4 
8.2 
1.1 


lo-'' 

10-* 
10-* 
10-* 
10-* 


En>0 
Eh>0 
Eh>0 
Eh>0 
0.13 


2.23 
2.61 
0.99 
0.18 

2.87 


3.34 
3.61 
1.83 
0.77 
3.77 



"For negative ions, aug-cc-pVXZ basis sets were used for the or- 
bital expansion. 

Table I: The total energy difference AE of CLDA from plain 
LDA the total negative effective repulsive charge Q„ (in e) 
and the IPs for selected atoms, molecules (top) and negative 
ions (bottom). IPs and EAs were calculated as the negative 
of the one-electron energies corresponding to the HOMO of 
the neutral system and the negative ion respectively. Basis 
set notation is explained in the caption of Fig. [S] For the 
neutral systems we compare with experimental values of the 
IP, while for the negative ions with experimental values of the 
EA of the corresponding neutral system. All energies are in 
eV. 



rect asymptotics are also obtained for larger systems like 
CO molecule as shown in Fig. [3] 

In the top of Table HI we show the IPs calculated with 
CLDA and with LDA, as the negative of the one-electron 
energy of the HOMO, i^n, for various atomic and molec- 
ular systems. In the inset of Fig. |3] we show that the 
IP value for CO (calculated in that way) is essentially 
independent of 9. 

The IPs from CLDA arc on average roughly 10% un- 
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Figure 3: The XC part, uxc, of the effective potential for CO 
molecule along the molecular axis. Green dotted lines indicate 
the correct ±l/r asymptotic behavior with r measured from 
the molecular center. In the inset, the En as a function of 
the SVD parameter 6 is shown with the horizontal blue line 
indicating the experimental value of the IP. 



derestimatcd. This divergence should be contrasted to 
the dramatic 40% errors of plain LDA. Given the se- 
vere underestimation of IPs and the fundamental gaps of 
solids by LDA and GGA, our approach offers a signifi- 
cant qualitative improvement. Contrary to LDA, nega- 
tive ions are predicted to be bound by CLDA as shown in 
the bottom of Table [B Even though EAs of neutral sys- 
tems (predicted as the negative of the one-electron energy 
of the HOMO of the corresponding negative ions), are 
underestimated compared to experiment by about 40%, 
it is nevertheless encouraging that qualitatively correct 
EAs can be obtained with CLDA. 

The differences of the total energies obtained with 
CLDA from those obtained with LDA are also shown 
in Table HI These differences are very small for all sys- 
tems, i.e., the additional constraint on the potential does 
not change the total energy. The total negative effective 
repulsive charge, (5„, is shown in Table H] as a measure of 
how well the positivity condition is fulfilled by the opti- 
mal potential. Although we did not manage to eliminate 
it completely, Qn is very small compared to the total ef- 
fective repulsive charge and does not affect the quality of 
our effective potential. 

As already mentioned, for CLDA we could have used 
the more sophisticated discontinuity correction of Ref. |36| 
instead of the SVD. To assess if this is required, we ap- 
plied the modified OEP equation (40) of Ref. [H to the Ne 
atom, adopting the same basis sets, and A = 10~^ as in 
the case of EXX. The obtained potential is almost on top 
of the CLDA potential shown in Fig. [2] while the IP, at 
19.6 eV, is only slightly improved compared to the result 
in Table m Thus, the choice of a simple SVD seemed to us 
reasonable for this first demonstration of our approach. 



V. DISCUSSION 

Until now, the main errors stemming from Sis were 
not in the total energy but resulted from deficiencies of 
the local potential. The main advantage of splitting the 
XC energy Exc ~ Ex + Eq and treating the exchange 
exactly is the cancellation of Sis and the quality of the 
KS potential, however at a computational cost compared 
with LDA/GGA, and an ensuing complicated description 
of correlation through non-local orbital functionals [401 . 
In fact an appropriate non-local Eq [p] of cost no- higher 
than EXX is not available yet. 

Attempting to address the SI problem in DFAs for fi- 
nite systems, we noted the ambiguity in the quantifica- 
tion of the SI error in the Hartree and XC energies. Still, 
it was possible to quantify the SI error in the KS poten- 
tial in a way that is unambiguous and independent of the 
error in the energy. Following our definition of the SI er- 
ror in the KS potential, we proposed two constraints for 
the repulsive part of the KS potential: the corresponding 
effective repulsive density must integrate to A'^ — 1 elec- 
trons and it must be a non-negative function everywhere. 
The constraint for the norm of the effective repulsive den- 
sity is a necessary property satisfied, for example, by the 
exact KS and by the EXX potentials but it is not a suf- 
ficient condition on its own. The norm together with the 
positivity constraint are sufficient conditions for the ab- 
sence of Sis from the potential, although the positivity 
constraint is probably too strong and not strictly satisfied 
by the exact KS potential and EXX. Since our treatment 
is still approximate, the ambiguity in dealing with Sis 
remains. 

Nevertheless, imposing these constraints constitutes a 
powerful method because the constraints restrict consid- 
erably the variational freedom of the effective potential 
while at the same time allowing for the accurate descrip- 
tion of EXX. When our constraints are applied to the 
EXX energy functional they produce potentials that are 
very close to the EXX potential obtained without any 
constraint. In particular, the constrained EXX potential 
in Fig.[T]preserves the shell structure of the numerical po- 
tential, by showing the corresponding bump at almost the 
same position. This excellent agreement, given the finite 
nature of the basis sets, shows that the constraints intro- 
duced here do not wash away the atomic shell structure 
from the optimal potential which is a subtle but essen- 
tial feature. The fact that this structure is absent from 
CLDA should be traced back to the approximate nature 
of the LDA functional. 

Finally, our constraints keep the scaling of computa- 
tional cost at the level of the corresponding DFA and 
eliminate the effects of Sis from the potential (where 
it matters) with a minimal increase of the total energy. 
For LDA, the constrained KS potentials have the correct 
asymptotic behavior and give significantly improved IPs 
over the unconstrained LDA results. At the same time, 
the description of XC has been kept together which has 
advantages, for example it exploits the cancellation of 



errors in Exc ^nd provides an improved description of 
electron-pair bonds |39j . 

It is accepted that the many-body SI error leads to spu- 
rious fractionally-charged atoms in the dissociation of a 
hetero-nuclear molecule. To prevent this spurious effect, 
the HXC potential must develop a barrier or a step be- 
tween the atoms, as shown by Perdew, in Ref. |4l|. These 
steps or barriers are related with the shell-structure in 
atoms; as we have seen in Figs. [1] [2l the latter is re- 
produced by our constraints in the case of EXX but not 
in LDA. Further investigation is necessary to examine 
whether our method leads to fractionally charged atoms 
in the dissociation of molecules. It is almost certain that 
this issue depends on the overall quality of the approxi- 
mate energy functional, and does not rely merely on the 
absence of Sis from the potential. 



Other interesting questions regarding our approach are 
its size consistency and how it applies to extended sys- 
tems where N and iV — 1 arc the same. All these ques- 
tions and challenges are related with the energetic cost 
to localize the XC density. In order to address them, our 
method needs to be extended to include a (SI) correc- 
tion energy term based on the correction of the effective 
potential. This is the focus of work which is in progress. 
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